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Abstract 



This paper concerns the reliability of a pair of coupled oscillators in response to fluctuating 
inputs. Reliability means that an input elicits essentially identical responses upon repeated pre- 
sentations regardless of the network's initial condition. Our main result is that both reliable and 
unreliable behaviors occur in this network for broad ranges of coupling strengths, even though 
individual oscillators are always reliable when uncoupled. A new finding is that at low input 
amplitudes, the system is highly susceptible to unreliable responses when the feedforward and 
feedback couplings are roughly comparable. A geometric explanation based on shear-induced 
chaos at the onset of phase-locking is proposed. 
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Introduction 

This paper, together with its companion paper Reliability of Coupled Oscillators II, contain a 
mathematical treatment of the question of reliability. Reliability here refers to whether a sys- 
tem produces identical responses when it is repeatedly presented with the same stimulus. Such 
questions are relevant to signal processing in biological and engineered systems. Consider, for 
example, a network of inter-connected neurons with some background activity. An external stim- 
ulus in the form of a time-dependent signal is applied to this neural circuitry, which processes 
the signal and produces a response in the form of voltage spikes. We say the system is reliable 
if, independent of its state at the time of presentation, the same stimulus elicits essentially identi- 
cal responses following an initial period of adjustment. That is, the response to a given signal is 
reproducible |l3l[5l|27llMl471[39l[M[30lEa[IlEa[3l. 

This study is carried out in the context of (heterogeneous) networks of interconnected oscilla- 
tors. We assume the input signal is received by some components of the network and relayed to 
others, possibly in the presence of feedback connections. Our motivation is a desire to understand 
the connection between a network's reliability properties and its architecture, i.e. its "circuit dia- 
gram," the strengths of various connections, etc. This problem is quite different from the simpler 
and much studied situation of uncoupled oscillators driven by a common input. In the latter, the 
concept of reliability coincides with synchronization, while in coupled systems, internal synchro- 
nization among constituent components is not a condition for reliability. 

To simplify the analysis, we assume the constituent oscillators are phase oscillators or circle 
rotators, and that they are driven by a fluctuating input which, for simplicity, we take to be white 
noise. Under these conditions, systems consisting of a single, isolated phase oscillator have been 
explored extensively and shown to be reliable; see, e.g., 1391 1341 . Our results are presented in a 
two-part series: 

Paper I. Two-oscillator systems 

Paper II. Larger networks 

The present paper contains an in-depth analysis of a 2-oscillator system in which the stimulus is 
received by one of the oscillators. Our results show that such a system can support both reliable 
and unreliable dynamics depending on coupling constants; a more detailed discussion of the main 
points of this paper is given later in the Introduction. Paper II extends some of the ideas developed 
in this paper to larger networks that can be decomposed into subsystems or "modules" so that the 
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inter-module connections form an acyclic graph. At the level of inter-module connections, such 
"acyclic quotient networks" have essentially feedforward dynamics; they are reliable if all the 
component modules are also reliable. Acyclic quotient networks are allowed to contain unreliable 
modules, however, and the simplest example of such a module is, as shown in this paper, an 
oscillator pair with nontrivial feedforward-feedback connections. In Paper II, we also explore 
issues surrounding the creation and propagation of unreliability in larger networks. 

With regard to aims and methodology, we seek to identify and explain relevant phenomena, and 
to make contact with rigorous mathematics whenever we can, hoping in the long run to help bring 
dynamical systems theory closer to biologically relevant systems. The notion of reliability, for 
example, is prevalent in neuroscience. With all the idealizations and simplifications at this stage of 
our modeling, we do not expect our results to be directly applicable to real-life situations, but hope 
that on the phenomenological level they will shed light on systems that share some characteristics 
with the oscillator networks studied here. Rigorous mathematical results are presented in the form 
of "Theorems," "Propositions," etc., and simulations are used abundantly. Our main results are 
qualitative. They are a combination of rigorous statements and heuristic explanations supported by 
numerical simulations and/or theoretical understanding. 

Content of present paper 

A motivation for this work is the following naive (and partly rhetorical) question: Are networks of 
coupled phase oscillators reliable, and if not, how large must a network be to exhibit unreliable 
behavior? Our answer to this question is that unreliable behavior occurs already in the 2-oscillator 
configuration in Diagram ([T]). Our results demonstrate clearly that such a system can be reliable 
or unreliable, and that both types of behaviors are quite prominent, depending in a generally pre- 
dictable way on the nature of the feedforward and feedback connections. Furthermore, we identify 
geometric mechanisms responsible for these behaviors. 

Referring the reader again to Diagram ([T]), three of the system's parameters are and a^, the 
feedforward and feedback coupling constants, and e, the amplitude of the stimulus. The following 
is a preview of the main points of this paper: 

1 . Lyapunov exponents as a measure of reliability. Viewing the stimulus-driven system as 
described by a stochastic differential equation and leveraging existing results from random 
dynamical systems theory, we explain in Sect. 2 why the top Lyapunov exponent (Amax) of 
the associated stochastic flow is a reasonable measure of reliability. In this interpretation, 
reliability is equated with Amax < 0, which is known to be equivalent to having random 
sinks in the dynamics, while unreliability is equated with Amax > 0, which is equivalent to 
the presence of random strange attractors. 

2. Geometry and zero-input dynamics. In pure feedforward and feedback configurations, i.e., 
at ttfb = or Off = 0, we identify geometric structures that prohibit unreliability. Our main 
result about zero-input systems is on phase-locking: in Sect. 3.3, we prove that for all 

in a broad range, the system is 1 : 1 phase-locked for either > or afb < ag depending 
on the relative intrinsic frequencies of the two oscillators. This phenomenon has important 
consequences for reliability. 
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3. Shear-induced chaos as main cause for unreliability at low drive amplitudes. Recent ad- 
vances in dynamical systems theory have identified a mechanism for producing chaos via 
the interaction of a forcing with the underlying shear in a system. The dynamical envi- 
ronment near the onset of phase-locking is particularly susceptible to this mechanism. We 
verify in Sect. 4.2 that the required "shearing" is indeed present in our 2-oscillator system. 
Applying the cited theory, we are able to predict the reliability or lack thereof for coupling 
parameters near the onset of phase-locking. At low drive amplitudes, this is the primary 
cause for unreliability. 

4. Reliability profile as a function ofag, Ofb and e. Via numerical simulations and theoretical 
reasoning, we deduce a rough reliability profile for the 2-oscillator system as a function 
of the three parameters above. With the increase in |afr|, |afb| and/or e, both reliable and 
unreliable regions grow in size and become more robust, meaning Amax is farther away from 
zero. The main findings are summarized in Sect. 4.3. 



1 Model and Formulation 
1.1 Description of the model 

We consider in this paper a small network consisting of two coupled phase oscillators forced by an 
external stimulus as shown: 

We assume the coupling is via smooth pulsatile interactions as in [l37ll45l [8l[T2ll. and the equations 
defining the system are given by 

^1 = ui + a^z{ei)g{e2)+ez{9{)l{t) , 
62 = UJ2 + as z{e2) g{9i) . 

The state of each oscillator is described by a phase, i.e., an angular variable 6'^ G §^ = M/Z, i = 
1,2. The constants ui and uj2 are the cells' intrinsic frequencies. We allow these frequencies to 
vary, but assume for definiteness that they are all within about 10% of 1. The external stimulus 
is denoted by el{t). Here I{t) is taken to be white noise, and e is the signal's amplitude. Notice 
that the signal is received only by cell 1. The coupling is via a "bump function" g which vanishes 
outside of a small interval (—6, b). On (—6, b), g is smooth, it is > 0, has a single peak, and satisfies 
J^^ g{6) dO = 1. The meaning of g is as follows: We say the ith oscillator "spikes" when 9i{t) = 0. 
Around the time that an oscillator spikes, it emits a pulse which modifies the other oscillator. The 
sign and magnitude of the feedforward coupling is given by ("forward" refers to the direction of 
stimulus propagation): > (resp. as < 0) means oscillator 1 excites (resp. inhibits) oscillator 

2 when it spikes. The feedback coupling, afb, plays the complementary role. In this paper, b is 
taken to be about ^, and as and afb are taken to be 0{1). Finally, z(9), often called the phase 
response curve [|44l [151 1811411231. measures the variable sensitivity of an oscillator to coupling and 
stimulus input. In this paper, we take z{6) = 2^(1 — cos(27r^^)) (see below). 
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Figure 1: Raster plot for an isolated oscillator. In each experiment, 75 trials are performed, and dots are 
placed at spike times. Nearly identical spike times are observed after a transient, indicating reliability. 

1.2 Neuroscience interpretations 

Coupled phase oscillators arise in many settings ll32l[35ll23l[T9ll45l[33l . Here, we briefly discuss 
their use in mathematical neuroscience. 

We think of phase oscillators as paradigms for systems with rhythmic behavior. Such models 
are often derived as limiting cases of oscillator models in two or more dimensions. In particular, 
the specific form of z(-) chosen here corresponds to the normal form for oscillators near saddle- 
node bifurcations on their limit cycles [8J. This situation is typical in neuroscience, where neural 
models with z(9) ^ 1 — cos(-) are referred to as "Type I." The pulse- or spike-based coupling 
implemented hy g{-) may also be motivated by the synaptic impulses sent between neurons after 
they fire action potentials (although this is not the only setting in which pulsatile coupling arises 
ESI [1 [m 1371 [361 [m |29l [311 [111). 

The general conclusions that we will present do not depend on the specific choices of and 
g{-), but rather on their qualitative features. Specifically, we have checked that our main results 
about reliability and phase locking are essentially unchanged when the z{-) function becomes 
asymmetric and the location of the g{-) impulse is somewhat shifted, as would correspond more 
closely to neuroscience Therefore, the behavior we find here can be expected to be fairly 
prototypical for pairs of pulse-coupled Type I oscillators. 

A standard way to investigate the reliability of a system of spiking oscillators — both in the 
laboratory and in simulations — is to repeat the experiment using a different initial condition 
each time but driving the system with the same stimulus el{t), and to record spike times in raster 
plots. Fig. [H shows such a raster plot for an isolated oscillator of the present type, as studied 
by [331 [HI. Note that, for each repetition, the oscillator produces essentially identical spike times 
after a transient, demonstrating its reliability. 

2 Measuring Reliability 

All of the ideas discussed in this section are general and are easily adapted to larger networks. 



5 



2.1 A working definition of reliability 



We attempt to give a formal definition of reliability in a general setting. Consider a dynamical 
system on a domain M (such as a manifold or a subset of Euclidean space). A signal in the form 
of bit) E M", t E [0, oo), is presented to the system. The response F(t) of the system to this signal 
is given by F{t) = F(xo, t, {i(s)}o<s<t). That is to say, the response at time t may, in principle, 
depend on xq E M, the initial state of the system when the signal is presented, and the values of 
the signal up to time t. 

In the model as described in Sect. 1.1, F(t) can be thought of as the pair {9i(t) , 92{t)) or 
\E'(6'i(t), 02{t)), the value of an observable at time t. 

We propose now one way to define reliability. Given a dynamical system, a class of signals X, 
and a response function F, we say the system is reliable if for almost all l E X and xo,Xq E M, 



Here || • || is a norm on the range space of F. We do not claim that this is the only way to capture 
the idea of reliability, but will use it as our operational definition. 

We point out some of the pitfalls of this definition: In practice, signals are never presented for 
infinite times, and in some situations, responses can be regarded as reliable only if the convergence 
above is rapid. By the same token, not all initial conditions are equally likely, leaving room for 
probabilistic interpretations. 

Finally, one should not expect unreliable responses to be fully random. On the contrary, as we 
will show in Sect. 2.2, they tend to possess a great deal of structure, forming what are known as 
random strange attractors. 

2.2 Reliability, Lyapunov exponents, and random attractors 

We discuss here some mathematical tools that can be used to quantify how reliable or unreliable a 
driven system is. With I{t) taken to be realizations of white noise, ^ can be put into the framework 
of a random dynamical system (RDS). We begin by reviewing some relevant mathematics yl|2l|. 
Consider a general stochastic differential equation 



In this general setting, Xt E M where M is a compact Riemannian manifold, the are indepen- 
dent standard Brownian motions, and the equation is of Stratonovich type. Clearly, Eq. ([21) is a 
special case of Eq. dS]): Xt describes the phases of the 2 oscillators, M = = x §^ (thus we 
have the choice between the Ito or Stratonovich integral), and k = 1. 

In general, one fixes an initial xq, and looks at the distribution of xt for t > 0. Under fairly 
general conditions, these distributions converge to a unique stationary measure fx, the density of 
which is given by the Fokker-Planck equation. Since reliability is about a system's reaction to 
a single stimulus, i.e., a single realization of the W^, at a time, and concerns the simultaneous 
evolution of all or large ensembles of initial conditions, of relevance to us are not the distributions 



\\F{xo,t,{t{s)}o<s<t) - F{x'Q,t,{L{s)}o<s<t)\\ as t ^ oo . 



k 




(3) 



1=1 
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of Xt hut flow-maps of the form Ft^ t^.^. Here ti < t2 are two points in time, lu is a sample 
Brownian path, and Ft^ t^.^{xt-^) = Xt^ where Xt is the solution of ([3]) corresponding to uj. A 
well known theorem states that such stochastic flows of dijfeomorphisms are well defined if the 
functions a(x) and b(x) in Eq. ^ are sufficiently smooth. More precisely, the maps Ft^^t2;u} 
are well defined for almost every uj, and they are invertible, smooth transformations with smooth 
inverses. Moreover, Ft^^^t^.^ and Fp^^^.^ are independent for ti < t2 < ^3 < ^4- These results allow 
us to treat the evolution of systems described by ([3]) as compositions of random, Ltd., smooth 
maps. Many of the techniques for analyzing smooth deterministic systems have been extended to 
this random setting. We will refer to the resulting body of work as RDS theory. 

Similarly, the stationary measure fx, which gives the steady-state distribution averaged over all 
realizations lu, does not describe what we see when studying a system's reliability. Of relevance 
are the sample measures {fXui}, which are the conditional measures of /i given the past. Here we 
think of to as defined for all t G (— oo, oo) and not just for t > 0. Then fi^ describes what one 
sees at t = given that the system has experienced the input defined by cu for all t < 0. This is 
an idealization: in reality, the input is presented on a time interval [—to, 0) for some to > 0. Two 
useful facts about these sample measures are 

(a) {F^t,o;uj)*fJ' — /Uo; as t — > oo, where (-F-t,o;(^)*Ai is the measure obtained by transporting fx 
forward by F^t,o;uj, and 

(b) the family {/x^} is invariant in the sense that (-fo,t;Lj)*(/U<^) = fJ-atiuj) where at{uj) is the time- 
shift of the sample path lu by t. 

Thus, if our initial distribution is given by a probability density p and we apply the stimulus corre- 
sponding to LU, then the distribution at time t is (-fo,i;(^)*P- For t sufficiently large, one expects in 
most situations that {FQ^t;u))*P is very close to (-fo,t;<^)*/^5 which by (a) above is essentially given by 
l^atiuj) ■ The time-shift by t of is necessary because by definition, fi^ is the conditional distribution 
of fi at time 0. 

Fig. [21 shows some snapshots of {FQ t;uj)*P for the coupled oscillator pair in the system (1) 
for two different sets of parameters. In the panels corresponding to t ^ 1, the distributions 
approximate fiat{u}) ■ In these simulations, the initial distribution p is the stationary density of Eq. ^ 
with e ~ 0.01, the interpretation being that the system is intrinsically noisy even in the absence 
of external stimuli. Observe that these pictures evolve with time, and as t increases they acquire 
similar qualitative properties depending on the underlying system. This is in agreement with RDS 
theory, which tells us in fact that the fiat(ui) obey a statistical law for almost all uj. Observe also the 
strikingly different behaviors in the top and bottom panels. We will follow up on this observation 
presently. First, we recall two mathematical results that describe the dichotomy. 

In deterministic dynamics, Lyapunov exponents measure the exponential rates of separation 
of nearby trajectories. Let Xmiixix) denote the largest Lyapunov exponent along the trajectory 
starting from x. Then a positive Amax for a large set of initial conditions is generally thought of 
as synonymous with chaos, while the presence of stable equilibria is characterized by A^ax < 0. 
For smooth random dynamical systems, Lyapunov exponents are also known to be well defined; 
moreover, they are nonrandom, i.e. they do not depend on uu. If, in addition, the invariant measure 
is erogdic, then Amax is constant almost everywhere in the phase space. In Theorem 1 below, we 
present two results from RDS theory that together suggest that the sign of A^ax is a good criterion 
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t = 20 



t = 50 



t = 500 



t = 1900 




(a) Random fixed point (Amax < 0) 
t = 20 t = 50 t = 500 t = 1900 




(b) Random strange attractor (Amax > 0) 

Figure 2: Snapshots of sample measures for Eq. ^ at various times in response to a single realization of the 
stimulus. Two different sets of parameters are used in (a) and (b). In (a), the sample measures converge to a 
random fixed point. In (b), the sample measures converge to a random strange attractor See Theorem 1. 



for distinguishing between reliable and unreliable behavior: 

Theorem 1. In the setting ofEq. (|21), let /i be an ergodic stationary measure. 

(1) (Random sinks) [|20ll If Xmax < 0, then with probability 1, fi^ is supported on a finite set of 
points. 

(2) (Random strange attractors) 1*241 IffJ- hcis a density and Amax > 0, then with probability 1, 
ficu is a random SRB measure. 

The conclusion of Part (1) is interpreted as follows: The scenario in which the support of 
/i^ consists of a single point corresponds exactly to reliability for almost every to as defined in 
Sect. 12.11 This was noted in, e.g., [|34l [30l . For the 2-oscillator system defined by Eq. the 
collapse of all initial conditions to a point is illustrated in Fig. Oa); notice that iXf^^^ continues to 
evolve as t increases. Even though in general //^ can be supported on more than one point when 
Amax < 0, this seems seldom to be the case except in the presence of symmetries. We do not know 
of a mathematical result to support this empirical observation, howeverQ In the rest of this paper, 

'Analysis of single-neuron recordings have revealed firing patterns which suggest the possible presence of random 
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we will take the liberty to equate Amax < with reliability. 

The conclusion of Part (2) requires clarification: In deterministic dynamical systems theory, 
SRB measures are natural invariant measures that describe the asymptotic dynamics of chaotic 
dissipative systems (in the same way that Liouville measures are the natural invariant measures for 
Hamiltonian systems). SRB measures are typically singular. They are concentrated on unstable 
manifolds, which are families of curves, surfaces etc. that wind around in a complicated way in 
the phase space [7]. Part (2) of Theorem 1 generalizes these ideas to random dynamical systems. 
Here, random (meaning cj-dependent) SRB measures live on random unstable manifolds, which 
are complicated families of curves, surfaces, etc. that evolve with time. In particular, in a system 
with random SRB measures, different initial conditions lead to very different outcomes at time t 
when acted on by the same stimulus; this is true for all t > 0, however large. It is, therefore, 
natural to regard Amax > and the distinctive geometry of random SRB measures as a signature of 
unreliability. 

In the special case where the phase space is a circle, such as in the case of a single oscillator, 
that the Lyapunov exponent A is < is an immediate consequence of Jensen's Inequality. In more 
detail, 

A(x) = lim -logF^^ (x) 

i-^oo t 

for typical uj by definition. Integrating over initial conditions x, we obtain 

A = / lim ilogF^^ (x) rfx = lim i / XogF'^ ^ {x) dx . 

The exchange of integral and limit is permissible because the required integrability conditions are 
satisfied in stochastic flows ETI . Jensen's Inequality then gives 

/ \ogF^^^Jx)dx < log / F^,.Jx)dx = 0. (4) 

The equality above follows from the fact that Fq is a circle diffeomorphism. Since the gap in 
the inequality in dH) is larger when Fq ^.^ is farther from being a constant function, we see that 
A < corresponds to F^^.^ becoming "exponentially uneven" as t ^ oo. This is consistent with 
the formation of random sinks. 

The following results from general RDS theory shed light on the situation when the system is 
multi-dimensional: 

Proposition 2.1. (see e.g. fll]) In the setting ofEq. assume /i has a density, and let - ■ ■ , A^} 

be the Lyapunov exponents of the system counted with multiplicity. Then 

(i) E^ Ai<0; 

(ii) J2i -^4 = if and only if Fs^t,uj preserves fifor almost all uj and all s < t; 
(Hi) i/^j Aj < 0, and Aj 7^ for all i, then /i^ is singular 

sinks with > 1 point [JJJ . 
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A formula giving the dimension of fi^ is proved in [|24ll under mild additional conditions. 

The reliability of a single oscillator, i.e. that A < 0, is easily deduced from Prop. 12. U ji has 
a density because the transition probabilities have densities, and no measure is preserved by all 
the Fs^t,uj because different stimuli distort the phase space differently. Prop. l2.lh ') and (ii) together 
imply that A < 0. See also, e.g., [ISOllMlllllMl- 

The remarks above concerning ji apply also to the 2-oscillator model in Eq. Q. (That fi has a 
density is explained in more detail in Sect. 3.2.) Therefore, Prop. |2.1t i) and (ii) together imply that 
Ai + A2 < 0. Here Ai = A^ax can be positive, zero, or negative. If it is > 0, then it will follow from 
Prop. I2.lt i) that A2 < 0, and by Prop. l2.lf iii). the ji^ are singular. From the geometry of random 
SRB measures, we conclude that different initial conditions are attracted to lower dimensional 
sets that depend on stimulus history. Thus even in unreliable dynamics, the responses are highly 
structured and far from uniformly distributed, as illustrated in Fig.|2b). 

Finally, we observe that since A^ax is nonrandom, the reliability of a system is independent of 
the realization uj once the stimulus amplitude e is fixed. 

3 Coupling geometry and zero-input dynamics 

3.1 Preliminary observations 

First we describe the flow ipt on the 2-torus defined by Eq. ^ when the stimulus is absent, i.e., 
when £ = 0. We begin with the case where the two oscillators are uncoupled, i.e. Ofj = afb = 0. In 
this special case, f^t is a linear flow; depending on whether uji/uj2 is rational, it is either periodic or 
quasiperiodic. Adding coupling distorts flow lines inside the two strips {|6'i| < h] and {\62\ < h}. 
These two strips correspond to the regions where one of the oscillators "spikes," transmitting a 
coupling impulse to the other (see Sect. 1.1). For example, if ag > 0, then an orbit entering 
the vertical strip containing = will bend upward. Because of the variable sensitivity of the 
receiving oscillator, the amount of bending depends on ag as well as the value of 6*2, with the effect 
being maximal near 6*2 = | and negligible near 6*2 = due to our choice of the function z{9). The 
flow remains linear outside of these two strips. See Fig.[3l 

Because the phase space is a torus and Eq. ^ does not have fixed points for the parameters 
of interest, ipt has a global section, e.g., Sq = {62 = 0}, with a return map Tq : Sq ^ So- The 
large-time dynamics of Lpt are therefore described by iterating Tq, a circle diffeomorphism. From 
standard theory, we know that circle maps have Lyapunov exponents < for almost all initial 
conditions with respect to Lebesgue measure. This in turn translates into Amax = for the flow Lpt- 

3.2 Reliability of two special configurations 

We consider next two special but instructive cases of Eq. namely the "pure feedforward" case 
corresponding to Ofb = and the "pure feedback" case corresponding to ag = 0. We will show 
that in these two special cases, the geometry of the system prohibits it from being unreliable: 

Proposition 3.1. For every £ > 0, 
(a) the system (|2]) has an ergodic stationary measure jj, with a density; 
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(b) (i) ifa^ = 0, then A^ax < 0; 
(ii) ifttfi = 0, then Amax = 0. 

We first discuss these results informally. Consider the case Ofb = 0. Notice that when e = 0, 
the time-t map of the flow generated by Eq. ^ sends vertical circles (meaning sets of the form 
{6i = c} where c is a constant) to vertical circles. As our stimulus acts purely in the horizontal 
direction and its magnitude is independent of 02, vertical circles continue to be preserved by the 
flow-maps Fs^t,L^ when e > 0. (One can also see this from the variational equations.) As is well 
known, A^ax > usually results from repeated stretching and folding of the phase space. Maps 
that preserve all vertical circles are far too rigid to allow this to happen. Thus we expect Amax < 
when Ofb = 0. In the pure feedback case = 0, the second oscillator rotates freely without input 
from either external stimuli or the first oscillator. Thus the system always has a zero Lyapunov 
exponent corresponding to the uniform motion in the direction of 02. 

Before proving Proposition 13. 1[ we recall some properties of Lyapunov exponents. Let Fq 
denote the stochastic flow on T^, and let be a stationary measure of the process. Recall that for 
a.e. uj, /x-a.e. a; e and every tangent vector v at x, the Lyapunov exponent 

\uj{x,v) = lim jlog |DFo,t;^(x)f I 

is well defined. Moreover, if Ai < A2 are the Lyapunov exponents at x, then 

lim ^log|det(DFo,t;^(x))| = Ai + A2 . (5) 

We say {vi, V2} is a Lyapunov basis at x if A<^(x, Vi) = Aj. It follows from Q that for such a pair 
of vectors, 

lim -\og\smZ{DFo^t.^{x)vi,DFo^t;uj{x)v2)\ = 0, (6) 

t— >oo t 
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that is, angles between images of vectors in a Lyapunov basis do not decrease exponentially fast. 
It is a standard fact that Lyapunov bases exist almost everywhere, and that any nonzero tangent 
vector can be part of such a basis. 

Proof of Proposition \3. 1\ First, we verify that Eq. Q has a stationary measure with a density. 
Because the white noise stimulus el{t) instantaneously spreads trajectories in the horizontal (6'i) 
direction, an invariant measure must have a density in this direction. At the same time, the de- 
terministic part of the flow carries this density forward through all of since flowlines make 
approximately 45 degree angles with the horizontal axis. Therefore the two-oscillator system has 
a 2-D invariant density whenever e > 0. 

We discuss the pure feedforward case. The case of = is similar and is left to the reader. 

As noted above, when Ofb = the flow FQ^t;u leaves invariant the family of vertical circles. 
This means that Fqi-^ factors onto a stochastic flow on theflrst circle. More precisely, if tt : = 

X — > is projection onto the first factor, then for almost every uj and every s < t, there is 
a diffeomorphism Fq^huj of with the property that for all a: G T^, Tc{Fo^t;ujix)) = Fo^t.^{Tc{x)). 
One checks easily that Fg^t,^ is in fact the stochastic flow corresponding to the system in which 
oscillator 2 is absent and oscillator 1 alone is driven by the stimulus. For simplicity of notation, for 
X e and tangent vector v at x, we denote 7r(x) by x and D7i(x) -vhy v. Furthermore, we let A^^ 
denote the Lyapunov exponent of Fs^t;uj. From Prop. 12. ll it follows that A^^ < almost everywhere. 

We now consider the Lyapunov exponents of system At fi-a.e. x G T^, we pick a vector v 
in the -direction. Since FQ^t;uj preserves vertical circles and /i has a density, an argument similar 
to that before Prop. 12.1 1 gives \u,{x, v) < 0. Let {u, v} be a Lyapunov basis at x. From Eq. Q and 
the relation between Fo^t;^, and Fs^t;u), we deduce that 



X^{x, u) = lim - log |DFo,t;^(x) ■ u\ = lim - log \DFo^t;oj{^) ■ u\ = A<^(x, u) < . 

□ 

It is likely that Amax is typically < in the case Ofb = 0. We have shown that one of the 
Lyapunov exponents, namely the one corresponding to the first oscillator alone, is < 0. The value 
of Amax therefore hinges on \uj{x, v) for vectors v without a component in the -direction. Even 
though this growth rate also involves compositions of random circle maps, the maps here are not 
i.i.d.: the kicks received by the second oscillator are in randomly timed pulses that cannot be put 
into the form of the white noise term in Eq. ([3]). We know of no mathematical result that guarantees 
a strictly negative Lyapunov exponent in this setting, but believe it is unlikely that Eq. ([3]) will have 
a robust zero Lyapunov exponent unless ag = 0. 

To summarize, we have shown that without recurrent connections, it is impossible for a two- 
oscillator system to exhibit an unreliable response. 

3.3 Phase locking in zero-input dynamics 

This subsection is about the dynamics of the coupled oscillator system when e = 0. Recall that the 
intrinsic frequencies of the two oscillators, Ui and uj2, are assumed to be ^ 1. Our main result is 
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Figure 4: The critical value as functions of . In both plots, the dashed line is the curve (see text), 
and the shaded regions are the parameters for which 1:1 phase-locking occurs. 



that if and 002 differ by a little, then in regions of the (ag, Ofb) -parameter plane in the form of a 
square centered at (0, 0), the two oscillators are 1:1 phase-locked for about half of the parameters. 
Two examples are shown in Fig. HI and more detail will be given as we provide a mathematical 
explanation for this phenomenon. (Phase locking of pairs of coupled phase oscillators is studied 
in e.g., |[T2l[T7l mulls EH. A primary difference is that we treat the problem on an open region of 
the (off, Ofb) -plane centered at (0, 0).) 

Let ipt be the flow on defined by ^ with e = 0. To study ipu we consider its return map 
T : Sfe — ^ Sfe where S;, = {62 = b}. Working with the section S;, (as opposed to e.g. Sq) simplifies 
the analysis as we will see, and substantively the results are not section dependent. Let p(T) be the 
rotation number of T. Since ui ^ ^012, it is natural to define p{T) in such a way that p(T) ^ 1 when 
ctff = ctfb = 0, so that p(T) may be interpreted as the average number of rotations made by the 
first oscillator for each rotation made by the second. It is well known that if p{T) is irrational, then 
ift is equivalent to a quasi-periodic flow via a continuous change of coordinates, while p(T) E Q 
corresponds to the presence of periodic orbits for cpf In particular, attracting fixed points of T 
correspond to attracting periodic orbits of cpt that are 1:1 phase-locked, "1:1" here referring to the 
fact that oscillator 2 goes around exactly once for each rotation made by oscillator 1 . 

We begin with the following elementary facts: 

Lemma 3.1. (a) With toi, uj2, and affixed and letting T = T^^^, the function ^ Ta^{6i) is 
strictly increasing for each 9i; it follows that p{Taj^^) is a nondecreasing functioi^of ai^. 
(b) With ui, as, and a^k fixed, 002 ^— > T^2{9i) is strictly decreasing for each 61, and p(T^2) ^ 
non-increasing function of 002. 

Analogous results hold as we vary and ui separately keeping the other three quantities fixed. 

Proof. This follows from the way each coupling term bends trajectories; see Sect. 3.1 and Fig.[3l 
We show (a); the rest are proved similarly. Consider two trajectories, both starting from the same 

-This is to allow for phase locking at rational values of p(Tafb)- 
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point on E;, but with different Ofb. They will remain identical until their ^2 -coordinates reach 1 — 6, 
as Ofb does not affect this part of the flow. Now at each point in the horizontal strip H = {1 — b < 
62 < 1 + 6}, the vector field corresponding to the larger Ofb has a larger horizontal component 
while the vertical components of the two vector fields are identical. It follows that the trajectory 
with the larger Ofb will be bent farther to the right as it crosses H . □ 

Our main result identifies regions in parameter space where T has a fixed point, corresponding 
to 1 : 1 phase-locking as discussed above. We begin with the following remarks and notation: (i) 
Observe that when uji = 102 and as = Ofb, we have T{x) = x for all x E Sf,; this is a consequence 
of the symmetry of z{9) about = \- (ii) We introduce the notation Acj = 1 — so that when 
flff = o,fb = 0, X — T{x) = Auj for all x, i.e., Ato measures the distance moved by the return map 
T under the linear flow when the two oscillators are uncoupled. Here, we have used T to denote 
not only the section map but its lift to M with T(l) = 1 in a harmless abuse of notation. 

We will state our results for a bounded range of parameters. For definiteness, we consider 
«fr, G [—2, 2] and .9 < cui < 1.1. The bounds for ag and Ofb are quite arbitrary. Once chosen, 
however, they will be fixed; in particular, the constants in our lemmas below will depend on these 
bounds. It is implicitly assumed that all parameters considered in Theorem 2 are in this admissible 
range. We do not view 6 as a parameter in the same way as tui, a;2, ag or Ofb, but instead of fixing 
it at we will take h as small as need be, and assume \g\ = 0{^) in the rigorous results to follow. 

Theorem 2. The following hold for all admissible {ui, uj2,as,a{h) and all b sufficiently small: 

(a) If UJ2 > uJi, then there exist = a^(co'i, 072, as) > as and £ = i{Au) > such that T has 
a fixed point for Ofb G [a^, + £] and no fixed point for Ofb < a^. 

(b) If UJ2 < uJi, then there exist = a^(ci;i, 0^2, as) < as and £ = £{Auj) > such that T has 
a fixed point for Ofb G [a^ — £, a^] and no fixed point for afb > a^. 

Moreover, \a^ — cls\ = O^Au); and for each Au 0, £ increases as b decreases. 

To prove this result, we need two lemmas, the proofs of both of which are given in the Ap- 
pendix. Define Aofb = — ag. 

Lemma 3.2. There exist bi and K > such that for all admissible {ui, uj2, as, ctfb), ifb < bi and 
Aofb < Kb-^Auo, then T(l - 6) < 1 - b. 

Lemma 3.3. There exist 62, C* > and xi G (0,1 — b) such that for all admissible {uJi,U2, as, apo), 
ifb < 62 cind Aflfb > CAuo, then T{xi) > xi. 

Proof of Theorem^ We prove (a); the proof of (b) is analogous. Let 071,^2 and as be given. 
Requirements on the size of b will emerge in the course of the proof. 

Observe first that with uj2 > uJi,T has no fixed point (and hence there is no 1:1 phase locking) 
when Ofb = Off. This is because T(x) = x when uj2 = coi and Ofb = as as noted earlier, and using 
Lemma [SlT b). we see that for 102 > <jJ\ and Ofb = ag, the graph of T is strictly below the diagonal. 

Keeping iO\,uj2 and as fixed, we now increase a^ starting from a^ = as. By Lemma [3Tt a). 
this causes the graph of T to shift up pointwise. As a^ is gradually increased, we let be the 
first parameter at which the graph of T intersects the diagonal, i.e. where there exists x* G Sb such 
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Figure 5: The rotation number p and Lyapunov exponent Amin as functions of Cfb- 



that T(x*) = X* — if such a parameter exists. Appealing once more to Lemma [3711 we see that 
p{T) < 1 for all flfb < a^, so T can have no fixed point for these values of Ofb. 

We show now that exists, and that the phase-locking persists on an interval of afb beyond 
a^. First, if b is small enough, then by Lemma [3^ T(l — b) < 1 — b for all afb < flff + Kb^'^Ato. 
Now if b is small enough that Kb^"^ > C where C is as in Lemma [331 then for afb G [aff + 
CAcu, aff + Kb~'^Aijj], T{xi) > Xi for some a^i < 1 — 6. For afb in this range, T maps the interval 
[xi, 1 — 6] into itself, guaranteeing a fixed point. It follows that (i) a^ exists and is < aff + CAu, 
and (ii) T has a fixed point for an interval of afb of length i > (Kb^^ — C)Auj. This completes the 
proof. □ 

As noted in the proof, for as long as uj2 7^ t^i, the lengths of the phase-locked intervals, i, can 
be made arbitrarily large by taking b small. On the other hand, if we fix b and shrink Auj, then 
these intervals will shrink. This is consistent with the phenomenon that the phase-locked region 
lies on opposite sides of the diagonal afb = aff when we decrease uj2 past uji, as shown in Fig.lH 

Instead of tracking the numerical constants in the proofs, we have checked numerically that for 
6 = ^, the pictures in Fig. |4] are quite typical, meaning about 50% of the parameters are phase 
locked. Specifically, for Auo up to about 10% and |aff | < 2, a^ — Off < 0.2, so that the a^-curve 
describing the onset of phase-locking is still quite close to the diagonal aff = afb. Also, for Au; 
as small as |%, the phase locked intervals have length £ > 4. These facts together imply that for 
parameters in the admissible range, the pictures are as shown in Fig.|4l 

Fig. [5] shows numerically-computed rotation numbers p and the rates of contraction to the 
corresponding limit cycle, i.e. the smaller Lyapunov exponent Amin of the flow ipf Notice that as 
afb increases past aj^, Amin decreases rapidly, so that the fixed point of T becomes very stable, a fact 
consistent with the large interval on which the system is 1 : 1 phase-locked. Furthermore, for the full 
range of |aff |, |afb| < 2 and 0.9 < oJi < 1.1, we find numerically that 0.53 < p{T) < 1.89. Phase- 
locking corresponding to rational p{T) with small denominators q {e.g., q = 3,4,5) is detected, 
but the intervals are very short and their lengths decrease rapidly with q. In other words, when 
the system is not 1:1 phase-locked - which occurs for about 50% of the parameters of interest - 
modulo fine details the system appears to be roughly quasi-periodic over not too large timescales. 
When the white-noise stimulus el{t) is added, these fine details will matter little. 
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Figure 6: The stretch-and-fold action of a kick followed by relaxation in the presence of shear. 

4 Reliable and Unreliable Behavior 

Numerical evidence is presented in Sect. 14.21 showing that unreliability can occur even when the 
stimulus amplitude is relatively small, and that its occurrence is closely connected with the onset of 
phase-locking in the zero-input system. A geometric explanation in terms of shear-induced chaos 
is proposed. Additionally, other results leading to a qualitative understanding of Amax as a function 
of Off, Ofb and e are discussed in Sect. 14.31 

4.1 A brief review of shear-induced chaos 

A rough idea of what we mean by "shear-induced chaos" is depicted in Fig.[6l An external force 
is transiently applied to a limit cycle (horizontal line), causing some phase points to move up and 
some to move down. Suppose the speeds with which points go around the limit cycle are height 
dependent. If the velocity gradient, which we refer to as "shear", is steep enough, then the bumps 
created by the forcing are exaggerated as the system relaxes, possibly causing the periodic orbit to 
fold. Such folding has the potential to lead to the formation of strange attractors. If the damping 
is large relative to the velocity gradient or the perturbation size, however, the bumps in Fig. [6] will 
dissipate before any significant stretching occurs. 

This subsection reviews a number of ideas surrounding the mechanism above. This mechanism 
is known to occur in many different dynamical settings. We have elected to introduce the ideas 
in the context of discrete-time kicking of limit cycles instead of the setting of Eq. (O because the 
geometry of discrete-time kicks is more transparent, and many of the results have been shown 
numerically to carry over with relatively minor modifications. Extensions to relevant settings are 
discussed later on in this subsection. A part of this body of work is supported by rigorous analysis. 
Specifically, theorems on shear-induced chaos for periodic kicks of limit cycles are proved in 
HOl HTl l42l l43l : it is from these articles that many of the ideas reviewed here have originated. 
Numerical studies extending the ideas in [ 4ni42ll to other types of underlying dynamics and forcing 
are carried out in [26J. For readers who wish to see a more in-depth (but not too technical) account 
of the material in this subsection, [26J is a good place to start. In particular. Sect. 1 in |[26l contains 
a fairly detailed exposition of the rigorous work in ll40ll4Tll42l[43l . 

Discrete-time kicks of limit cycles 

We consider a flow Lpt in any dimension, with a limit cycle 7. Let Tq < Ti < T2 < ■ ■ • be a 
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Figure 7: Geometry of folding in relation to the M^'"'-foIiation. Shown are a segment 70 C 7, its image after 
one kick, and two of the subsequent images under the flow. 

sequence of kick times, and let k„, n = 0, 1, 2, ■ ■ ■ , be a sequence of kick maps (for the moment 
Kn can be any transformation of the phase space). We consider a system kicked at time T„ by 
with ample time to relax between kicks, i.e., T„+i — T„ should not be too small on average. 

Central to the geometry of shear-induced chaos is the following dynamical structure of the 
unforced system: For each x G 7, the strong stable manifold W^^{x) of through x is defined to 
be 

W^x) = {y : lim d{ipt{x),ipt{y)) = 0} . 

These codimension 1 submanifolds are invariant under the flow, meaning ipt carries W^^{x) to 
W^^{'^t{'^)). In particular, if r is the period of the limit cycle, then lPt-{W^^{x)) = W^'^{x) for 
each X. Together these manifolds partition the basin of attraction of 7 into hypersurfaces, forming 
what is called the strong stable foliation. 

Fig. U\ shows a segment 70 C 7, its image 'Jq = /t(7) under a kick map k, and two images 
of 'Jq under Lfnr and ipmr for n > m G Z+. If we consider integer multiples of r, so that the 
flow-map carries each 14^ '"^ -leaf to itself, we may think of it as sliding points in 'Jq toward 7 along 
l^'^'^-leaves. (For t that are not integer multiples of r, the picture is similar but shifted along 7.) 
The stretching created in this combined kick- and- slide action depends both on the geometry of the 
H^'^^-foliation and on the action of the kick. Fig. |7] sheds light on the types of kicks that are likely 
to lead to folding: A forcing that drives points in a direction roughly parallel to the 1^** -leaves will 
not produce folding. Nor will kicks that essentially carry one W^^'^-leaf to another, because such 
kicks necessarily preserve the ordering of the leaves. What causes the stretching and folding is the 
variation along 7 in how far points a; G 7 are kicked as measured in the direction transverse to the 
W^^-leawes; we think of this as the "effective deformation" of the limit cycle 7 by the kick. 

To develop a quantitative understanding of the factors conducive to the production of chaos, it 
is illuminating to consider the following linear shear model [|46ll4TII : 

e = 1 + ay, 

y = -\y + A sin(27r^^) ■ ^^=0 ^ ' ^T) . ^ ^ 

Here, 6 E'B^ and ?/ G M are phase variables, and cr, A, A are constants^! We assume for definiteness 
that a and A are > 0, so that when A = 0, the system has a limit cycle 7 at {?/ = 0}. Its VT'^^-leaves 
are easily computed to be straight lines with slope = —K When A ^ the system is kicked 

^In this subsection, A denotes the damping constant in Eq. (Q. 
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periodically with period T. For this system, it has been shown that the shear ratio 
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is an excellent predictor of the dynamics of the system. Roughly speaking, if Amax denotes the 
largest observed Lyapunov exponent, then 

(a) if the shear ratio is sufficiently large, Amax is likely to be > 0; 

(b) if the shear ratio is very small, then Amax is slightly < or equal to 0; 

(c) as the shear ratio increases from small to large, Amax first becomes negative, then becomes 
quite irregular (taking on both positive and negative values), and is eventually dominated by 
positive values. 

To get an idea of why this should be the case, consider the composite kick-and-slide action in Fig.|7] 
in the context of Eq. (|7]). The time-T map of Eq. (|7]) is easily computed to be 



When the contraction in y is sufficiently strong, the first component of this map gives a good 
indication of what happens in the full 2-D system. As an approximation, define fxid) = 9^ and 
view fx as a map of 7 to itself. When the shear ratio is large and (1 — e^'^^) is not too small, 
\f^\ is quite large over much of 7, and the associated expansion has the potential to create the 
positive Amax mentioned in (a). At the opposite extreme, when the shear ratio is very small, fx 
is a perturbation of the identity; this is the scenario in (b). Interestingly, it is for intermediate 
shear ratios that fx tends to have sinks, resulting in Amax < for the 2-D system. The 1-D map 
/ = limT^oo /t is known variously as the phase resetting curve or the singular limit. It is used 
heavily in [|40l |4T1 l42l |43 1 to produce rigorous results for large T. 

We now return to the the more general setting of Lpt with discrete kicks k„, and try to interpret 
the results above as best we can. To make a meaningful comparison with the linear shear flow, 
we propose to first put our unforced system in "canonical coordinates," i.e., to reparametrize the 
periodic orbit 7 so it has unit speed, and to make the kick directions perpendicular to 7 — assuming 
the kicks have well defined directions. In these new coordinates, sizes of vertical deformations 
make sense, as do the idea of damping and shear, even though these quantities and the angles 
between iy**-manifolds and 7 all vary along 7. Time intervals between kicks may vary as well. 
The general geometric picture is thus a distorted version of the linear shear flow. We do not believe 
there is a simple formula to take the place of the shear ratio in this general setting; replacing the 
quantities a, A and A by their averages is not quite the right thing to do. We emphasize, however, 
that while system details affect the numerical values of Amax and the amount of shear needed to 
produce Amax > 0, the fact that the overall trends as described in (a)-(c) above are valid has been 
repeatedly demonstrated in simulation; see, e.g., [ 261 . 



Vt = 



6t 



9 + T + ^-[A sin(27re) + y] ■ (l - e"^^) (mod 1) , 
e"^^- {y + Asm{27Te)) . 



(9) 
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Generalization to stochastic forcing 

We now replace the discrete-time kicks in the discussion above by a directed (degenerate) 
continuous stochastic forcing, i.e. by a term of the form V{x)dWt where V is a vector field and 
dWt is white noise. By Trotter's product formula, the dynamics of the resulting stochastic flow can 
be approximated by a sequence of composite maps of the form (/^a* o KAt,Lu where At is a small 
time step, KAt,u; are kick maps (of random sizes) in the direction of V, and (^At is the unforced flow. 
Most of the time, the maps ^At^oj have negligible effects. This is especially the case if the size of 
V is not too large and damping is present. Once in a while, however, a large deviation occurs, 
producing an effect similar to that of the discrete-time kicks at times T„ described above. Cast in 
this light, we expect that the ideas above will continue to apply - albeit without the factors in the 
shear ratio being precisely defined. 

We mention two of the differences between stochastic forcing and periodic, discrete kicks. Not 
surprisingly, stochastic forcing gives simpler dependence on parameters: Amax varies smoothly, 
irregularities of the type in (c) above having been "averaged out". Overall trends such as those 
in (a)-(c) tend to be unambiguous and more easily detected than for deterministic kicks. Second, 
unlike periodic kicks, very small forcing amplitudes can elicit chaotic behavior without j being 
very large; this is attributed to the effects of large deviations. 

Generalization to quasi-periodic flows 

We have chosen to first introduce the ideas above in the context of limit cycles where the 
relevant geometric objects or quantities (such as a, A and W^'^) are more easily extracted. These 
ideas apply in fact to flows on a torus that are roughly "quasi-periodic" - meaning that orbits may 
or may not be periodic but if they are, the periods are large - provided the forcing, stochastic or 
discrete, has a well-defined direction as discussed earlier. The main difference between the quasi- 
periodic setting and that of a limit cycle is that ly^'^-leaves are generally not defined. A crucial 
observation made in [26] is that since folding occurs in finite time, what is relevant to the geometry 
of folding is not the usual iy**-foliation (which takes into consideration the dynamics as t ^ oo) 
hut finite -time strong stable foliations. Roughly speaking, a time-t W^^'^-leaf is a curve segment 
(or submanifold) that contracts in the first t units of time. For the ideas above to apply, we must 
verify that time-t iy*''-manifolds exist, have the characteristics of a large shear ratio, and that t is 
large enough for the folding to actually occur. If these conditions are met, then one can expect 
shear-induced chaos to be present for the same reasons as before. 

4.2 Phase-locking and unreliability in the two-cell model 

We now return to reliability questions of Eq. Q. In this subsection and the next, numerical data 
on Ainax as functions of as, afb and e are discussed. Our aim is to identify the prominent features 
in these numerical results, and to propose explanations for the phenomena observed. 

The present subsection is limited to phenomena related to the onset of phase-locking, which we 
have shown in Sect. 3.3 to occur at a curve in {as, afb)-space that runs roughly along the diagonal. 
We will refer to this curve as the a^-curve. Fig. [8] shows Amax as a function of as and Ofb for 
several choices of parameters, with a^-curves from Fig. |4] overlaid for ease of reference. In the top 
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e=0.2, co^=1 , (02=0.97 e=0.2, (o^=1 , ©2=1 .05 




(b) e = 0.8 

Figure 8: Maximum Lyapunov exponent Amax versus coupling strengths in the two-cell network. In all plots, 
we use uji = 1. The a^-curve from Fig.|4]is overlaid. All computed Amax shown here have standard errors of 
< 0.002 as estimated by the method of batched means. By the Central Limit Theorem, this means the actual 
Amax should lie within ^ 2.5 x 0.002 = 0.005 of the computed value with > 99% probability. Remark on 
plots: We have chosen the dynamic range in shading the figures to allow meaningful comparison of figures; 
a side effect is that some contour lines may not be visible. We always indicate the actual range of values 
through explicit labels. 
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0, 



Figure 9: Rotation of vectors in backwards time. Here, wi = 1, u;2 = 1.1, agf = 1, Ofb = 1.5. 

two panels, where e is very small, evidence of events connected with the onset of phase-locking is 
undeniable: definitively reliable (Amax < 0) and definitively unreliable (Amax > 0) regions are both 
present. Continuing to focus on neighborhoods of the a^-curves, we notice by comparing the top 
and bottom panels that for each (ag, afb), the tendency is to shift in the direction of unreliability 
as 6 is increased. We will argue in the paragraphs to follow that these observations are entirely 
consistent with predictions from Sect. 4.1. 

Shearing mechanisms 

For concreteness, we consider the case uj2 > uji, and consider first parameters at which the 
unforced system has a limit cycle, i.e., for each ag G [—1.5, 1.5], we consider values of Ofb that are 

> and not too far from a^. From Sect. 4.1, we learn that to determine the propensity of the 
system for shear-induced chaos, we need information on (i) the geometry of the limit cycle, (ii) the 
orientation of its ly^-manifolds in relation to the cycle, and (iii) the effective deformation due to 
the forcing. 

The answer to (i) is simple: As with all other trajectories, the limit cycle is linear with slope 

> 1 outside of the two corridors \9i\ < b and \92\ < b, where it is bent; see Fig.[3l As for (ii) and 
(iii), we already know what happens in two special cases, namely when ag = or Ofb = 0. As 
discussed in Sect. 3.2, when Ofb = 0, vertical circles are invariant under cpt. Since l^'^'^-leaves are 
the only manifolds that are invariant, that means the l^'^'^-manifolds are vertical. We noted also that 
the forcing preserves these manifolds. In the language of Sect. 4.1, this means the forcing creates 
no variation transversal to ly^'^-leaves: the ordering of points in this direction is preserved under 
the forcing. Hence shear-induced chaos is not possible here, and not likely for nearby parameters. 
A similar argument (which we leave to the reader) applies to the case = 0. From here on, we 
assume as, Ofb are both away from 0. 

We now turn to a treatment of (ii) when as, ctfb are both away from 0, and claim that ly^^'^-leaves 
generally have a roughly diagonal orientation, i.e., they point in a roughly southwest-northeast 
(SW-NE) direction. To find this orientation, we fix a point p on the limit cycle 7, and any nonzero 
tangent vector v at p (see Fig. |9l). We then flow backwards in time, letting p_t = f-tip) and 
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Figure 10: Strong stable directions along limit cycles. In (a), afb = + 0.1; in (b), Ofb = + 0.2. 
Additionally in (a), horizontal lines with arrows indicate the impact of the forcing at various points along 
the cycle; the variable impact is due to the z-function. 



v^t = D(p_t{p)v- The strong stable direction at p^t is the limiting direction of v^t-nr (t = period 
of 7) as n — > 00. Exact orientations of H^*'^ -leaves depend on as, afb and are easily computed 
numerically. 

The following simple analysis demonstrates how to deduce the general orientation of the W'^'^- 
leaves in the two-oscillator system from the signs of its couplings. For definiteness, we consider 
as > 0, so that is also positive and slightly larger than a^. Here, a typical situation is that 
if we identify the phase space with the square [0, 1] x [0, 1], then the limit cycle crosses the right 
edge {1} X [0, 1] in the bottom half, and the bottom edge [0, 1] x {0} in the right half (as shown 
in Fig. |9l). Let A, B and C be as shown, and consider a point p at A. Flowing backwards, suppose 
it takes time to reach point B, and time tc to reach point C . We discuss how V-t changes as we 
go from A to C . The rest of the time the flow is linear and V-t is unchanged. 
From A to B: Compare the backward orbits of p and p', where p' = p + kv and /c > is thought of 
as infinitesimally small. As these orbits reach the vertical strip {g > 0}, both are bent downwards 
due to as > 0. However, the orbit of p' is bent more because the function z(6) peaks at 6* = 1/2 
(see Fig.|9l). Thus, v^tg = v + {0, —5i) for some 5i > 0. 

From B to C: Continuing to flow backwards, we see by an analogous argument that as the two 
orbits cross the horizontal strip {g > 0}, both are bent to the left, and the orbit of p' is bent more. 
Therefore, V-tc = ""-ts + (— (^2, 0) for some ^2 > 0. 

Combining these two steps, we see that each time the orbit of p_t goes from Ato C, a vector 
of the form {—62,— Si) is added to v^f We conclude that as t 00, the direction of v^t is 
asymptotically SW-NE. Moreover, it is a little more W than S compared to the limit cycle because 
V-t must remain on the same side of the cycle. 

Numerical computations of strong stable directions are shown in Fig. \T0\ The plot in (a) is 
an example of the situation just discussed. The plot in (b) is for as, Ofb < 0, for which a similar 



22 



analysis can be carried out. Notice how small the angles are between the limit cycle and its W'^'^- 
leaves; this is true for all the parameter sets we have examined where Ofb is close to a^. Recall 
from Sect. 4.1 that it is, in fact, the angles in canonical coordinates that count. Since the limit 
cycle is roughly diagonal and the forcing is horizontal, putting the system in canonical coordinates 
will not change these angles by more than a moderate factor (except in one small region in picture 
(a)); i.e., the angles will remain small. 

Finally, we come to (iii), the deformation due to the forcing. Given that the forcing is in the 
horizontal direction and its amplitude depends on 9i (it is negligible when 9i ^ and has maximal 
effect when 9i ^ i), it causes "bump-like" perturbations transversal to the W^''* -manifolds (which 
are roughly SW-NE) with a geometry similar to that in Fig.|7l see Fig. fTOt a). 

This completes our discussion of the limit cycle case. We move now to the other side of the 
a^-curve, where the system is, for practical purposes, quasi-periodic (but not far from periodic). 
As discussed in the last part of Sect. 4.1, the ideas of shear-induced chaos continue to apply, 
with the role of iy**-leaves now played by finite-time stable manifolds. Since these manifolds 
change slowly with ag and afb, it can be expected - and we have checked - that they continue to 
make small angles with flowlines. Likewise, the forcing continues to deform flowlines by variable 
amounts as measured in distances transversal to finite-time stable manifolds. 

We conclude that when as, are both away from 0, the geometry is favorable for shear- 
induced stretching and folding. Exactly how large a forcing amplitude is needed to produce a 
positive Amax depends on system details. Such information cannot be deduced from the ideas 
reviewed in Sect. 4. 1 alone. 

Reliability-unreliability interpretations 

We now examine more closely Fig. [8l and attempt to explain the reliability properties of those 
systems whose couplings lie in a neighborhood of the a^-curve. The discussion below applies to 
|afj| > about 0.3. We have observed earlier that for ag or Ofb too close to 0, phase-space geometry 
prohibits unreliability. 

Consider first Fig.[8ta), where the stimulus amplitude e = 0.2 is very weak. Regions showing 
positive and negative Lyapunov exponents are clearly visible in both panels. Which side of the a^- 
curve corresponds to the phase-locked region is also readily recognizable to the trained eye (lower 
triangular region in the picture on the left and upper triangular region on the right; see Fig. 4). 

We first explore the phase-locked side of the a^-curve. Moving away from this curve, Amax 
first becomes definitively negative. This is consistent with the increased damping noted in Sect. 
3.2; see Fig. 6(b). As we move farther away from the a^-curve still, Amax increases and remains 
for a large region close to 0. Intuitively, this is due to the fact that for these parameters the limit 
cycle is very robust. The damping is so strong that the forcing cannot (usually) push any part of 
the limit cycle any distance away from it before it is brought back to its original position. That is 
to say, the perturbations are negligible. With regard to the theory in Sect. 4.1, assuming a remains 
roughly constant, that Amax should increase from negative to as we continue to move away from 

is consistent with increased damping; see (b) and (c) in the interpretation of the shear ratio. 

Moving now to the other side of the a^-curve, which is essentially quasi-periodic, regions of 
unreliability are clearly visible. These regions in fact begin slightly on the phase-locked side of 
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Figure 11: Folding action caused by white noise forcing and shear near the limit cycle (with afb > a^). At 
t = 0, the curve shown is the lift of the limit cycle 7 to M^. The remaining panels show lifts of the images 
Fo,t,uj{l) at increasing times. The parameters are wi = 1, 002 = 1-05, ag = 1, afb = 1.2, and £ = 0.8. Note 
that it is not difficult to find such a fold in simulations: very roughly, 1 out of 4 realizations of forcing gives 
such a sequence for t € [0, 5]. 



the curve, where a weakly attractive limit cycle is present. We have presented evidence to support 
our contention that this is due to shear-induced chaos, or folding of the phase space. The fact 
that Amax is more positive before the limit cycle is born than after can be attributed to the weaker- 
to-nonexistent damping before its birth. Thus the general progression of Amax from roughly 
to definitively negative to positive as we cross the a^-curve from the phase-locked side to the 
quasi-periodic side is altogether consistent with scenarios (a)-(c) in Sect. 4. 1 together with the 
observations in the paragraph on stochastic forcing. 

We point out that the unreliability seen in these panels is fairly delicate, perhaps even unex- 
pected a priori for the smaller values of ag and a^, such as 0.3 < |afj|, |afb| < 0.8: The bending 
of the flow lines is rather mild at these smaller coupling parameters. Moreover, we know that no 
chaotic behavior is possible at £ = 0, and the stimulus amplitude oi e = 0.2 in the top panels 
is quite small. Recall, however, that the stimulus is a fluctuating white noise, and e gives only 
an indication of its average amplitude. As noted in Sect. 4.1, we believe the unreliability seen is 
brought about by an interaction between the large fluctuations in the stimulus presented and the 
shearing in the underlying dynamics. 

In Fig. [Htb), the stimulus amplitude is increased to e = 0.8. A close examination of the plots 
shows that near to and on both sides of the a^-curve, Amax has increased for each parameter pair 
(ofj, Ofb), and that the reliable regions are pushed deeper into the phase-locked side compared to the 
top panels. This is consistent with the shear ratio increasing with forcing amplitude as predicted in 
Sect. 4.1. 

This completes our discussion in relation to Fig. [8l 

To complement the theoretical description of the geometry of folding given in Sect. 4.1, we 
believe it is instructive to see an actual instance of how such a fold is developed when the system 
Q is subjected to an arbitrary realization of white noise. A few snapshots of the time evolution of 
the limit cycle under such a forcing is shown in Fig.[TT] Notice that at the beginning, the combined 
action of the coupling and forcing causes the curve to wriggle left and right in an uncertain manner, 
but once a definitive kink is developed (such as at t = 3), it is stretched by the shear as predicted. 
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Figure 12: Amax as function of Ofb and e, with uji = I, uj2 = 1-1, and ag = 1. 

4.3 Further observations on parameter dependence 

We now proceed to other observations on the dependence of A^ax on network parameters. Our aim 
is to identify the salient features in the reliability profile of the oscillator system (|2l) as a function 
of as, Ofb and e, and to attempt to explain the phenomena observed. Our observations are based on 
plots of the type in Figs. [8] and [121 Some of the explanations we venture below are partial and/or 
speculative; they will be so indicated. 

1. Triple point: This phenomenon is the main topic of discussion in Sect. 4.2. Fig. [T2] shows a 
different view of the parameter dependence of Amax- At about Ofb = 1.4, which is near 
for the parameters used, both positive and negative Lyapunov exponents are clearly visible 
for very small £ in a manner that is consistent with Fig. [8] (even though the cu's differ slightly 
from that figure). We note again that this region, which we refer to as the "triple point," is 
an area of extreme sensitivity for the system, in the sense that the system may respond in a 
definitively reliable or definitively unreliable way to stimuli of very small amplitudes, with 
the reliability of its response depending sensitively on coupling parameters. 

2. Unreliability due to larger couplings: Fig. [8] and other plots not shown point to the occur- 
rence of unreliability when |afj| and |afb| are both relatively large. We are referring here 
specifically to the "off-diagonal" regions of unreliability (far from the a^-curve) in Fig. [8l 
This phenomenon may be partly responsible for the unreliability seen for the larger values 
of Off and afb on the diagonal as well; it is impossible to separate the effects of the different 
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mechanisms. 

We do not have an explanation for why one should expect Amax > for larger |aff| and |afb| 
aside from the obvious, namely that tangent vectors are more strongly rotated as they cross 
the strips {\9i\ < b}, making it potentially easier for folding to occur. But folding does not 
occur at £ = in spite of this larger bending. We believe the difference between the two 
situations is due to the following noise-assisted mechanism: For p G and a tangent vector 
u at p, let us say u is positively oriented with respect to flowlines if starting from the direction 
of the flow and rotating counter-clockwise, one reaches u before reaching the direction of 
the backward vector field. Without external forcing, if u is positively oriented, Dipt{p) ■ u 
will remain positively oriented for all t, because these vectors cannot cross flowlines. Now, 
in order for folding to occur, as in the formation of a horseshoe, the flow-maps must reverse 
the orientations of some tangent vectors. Even though larger values of \as\ and |afb| mean 
that tangent vectors are more strongly rotated, a complete reversal in direction cannot be 
accomplished without crossing flowlines. A small amount of noise makes this crossing 
feasible, opening the door (suddenly) to positive exponents. 

3. The ejfects of increasing e (up to around e = 2): 

(a) Unreliable regions grow larger, and Amax increases: A natural explanation here is that 
the stronger the stimulus, the greater its capacity to deform and fold the phase space - 
provided such folding is permitted by the underlying geometry. Because of the form 
of our stimulus, however, too large an amplitude simply pushes all phase points toward 
{6i = 0}. This will not lead to Amax > 0, a fact supported by numerics (not shown). 

(b) Reliable regions grow larger, and the responses become more reliable: As e increases, 
the reliable region includes all parameters (ag, afb) in a large wedge containing the 
Ofb = axis. Moreover, in this region Amax becomes significantly more negative as e 
increases. We propose the following heuristic explanation: In the case of a single oscil- 
lator, if we increase the amplitude of the stimulus, Amax becomes more negative. This 
is because larger distortions of phase space geometry lead to more uneven derivatives 
of the flow-maps Fs,t;a;, which in turn leads to a larger gap in Jensen's Inequality (see 
the discussion before Prop. 2.1). For two oscillators coupled as Eq. Q, increasing e 
has a similar stabilizing effect on oscillator 1 . Feedback kicks from oscillator 2 may 
destabilize it, as happens for certain parameters near the a^-curve. However, if 

is large enough and enhanced by a large e, it appears that the stabilizing effects will 
prevail. 
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5 Conclusions 



We have shown that a network of two pulse-coupled phase oscillators can exhibit both reliable and 
unreliable responses to a white-noise stimulus, depending on the signs and strengths of network 
connections and the stimulus amplitude. Specifically: 

1. (a) Dominantly feedforward networks are always reliable, and they become more reliable 

with increasing input amplitude. 

(b) Dominantly feedback networks are always neutral to weakly reliable. 

2. Whenfeedback and feedforward coupling strengths are comparable, whether both are nega- 
tive (mutually inhibitory) or positive (mutually excitatory), we have observed the following: 

(a) For weak input amplitudes, the system is extremely sensitive to coupling strengths, with 
substantially reliable and unreliable configurations occurring in close proximity. This 
phenomenon is explained by mechanisms of phase locking and shear-induced chaos. 

(b) For stronger input amplitudes, the system is typically unreliable. 

For weak stimuli, the most reliable configurations are, in fact, not dominantly feedforward but 
those with comparable feedforward-feedback couplings. 

We expect these results to be fairly prototypical for certain pulse-coupled phase oscillators, 
such as those with "Type I" phase response curves that frequently occur in neuroscience. Under- 
standing the effects of qualitatively different coupling types or oscillator models is an interesting 
topic for future study. 

Another natural extension is to larger networks. This is the topic of the companion paper [25^, 
where we use the present two-oscillator system as an embedded "module" to illustrate how unreli- 
able dynamics can be generated locally and propagated to other parts of a network. 
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Anne-Marie Oswald for helpful discussions over the course of this project. 



Here we prove the two lemmas needed for Theorem 2. 

Let c^i, Off , Atu and Aofb be given, with Aa;, Aofb > 0. To study the system where uj2 is defined 
by Ac<j = 1 — ^ and afb = ag + Aofb, we will seek to compare trajectories for systems with the 
following parameter sets: 



APPENDIX: Proof of Theorem 2 



System A : ag, Ofb 
System B : ag, Ofb 
System C : ag, Ofb 



= Off, CUi, U2=Ui 

= Off UJi^ UJ2 = UOi + Au! ■ UJ2 

= Off + AOfb, UJi, UJ2 = UJi + ACO ■ UJ2 
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Figure 13: Values used in proving Lemma ll!2l 



That is, System C is the system of interest, and Systems A and B are used to help analyze System 
C. We introduce also the following notation: H and V denote the horizontal and vertical strips of 
width 2b centered at integer values of 02 and 9i. We will work in instead of T^. 

Lemma \3.2i There exist hi and K > such that for all admissible {ui, uj2, a^, Ofb), ifb < hi and 
Aafb < Kb-^Auj, then T(l - 6) < 1 - b. 

Proof. For each of the 3 parameter sets above, two orbits are considered: Orbit 1 starts from 
(^1, 6*2) = (1 — b,b) E Sfe and runs forward in time until it meets = {62 = I — b}; orbit 2 
starts from (2 — 6, 1 + 6) G Si+t and runs backward in time until it meets Si_6. We need to prove 
that for System C, under the conditions in the lemma, the end point of orbit 1 lies to the left of the 
end point of orbit 2 (as shown in Fig. [13]). This is equivalent to r(l — 6) < 1 — 6. 

For System A, orbits 1 and 2 meet, since for this set of parameters, T(x) = x for all x as noted 
in Sect. 3.3. Comparing Systems A and B, since the vector field for System B has greater slope 
everywhere, and outside of if U ^ it has slope ^ > 1, we conclude that for System B the end 
point of orbit 1 lies to the left of the end point of orbit 2, with a separation h > Auj/2; see Fig.[T3l 

Next, we compare Systems B and C. Orbit 1 for the two systems is identical, since the equation 
outside of H does not involve Ofb. Orbit 2, however, differs for the two systems. To estimate by 
how much, we compare a, the distance from the end point of orbit 2 to 6^1 = 2 for System B, and 
a', the corresponding distance for System C as marked in Fig. [131 First, there exist bi and ki > 
such that for ^1 G (2 - 56i, 2), we have ^(^i) < A;i(2 - ^i)^ and |2'(^i)| < 2^1(2-^1). Shrinking 
bi further if necessary, we have that for b < bi, orbit 2 has slope > 1/2 everywhere and therefore 
the entire orbit, for both Systems B and C, lies within the region H D {2 — 5b < 9i < 2 — b}. The 
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next step is to apply Gronwall's Lemma to a system that incorporates both Systems B and C. Since 
92it) is identical for the two systems in the relevant region, we may write the equations as 

01 = -uji - {as + S)z{ei)g{t) 
6 = 0, 

where 6 = corresponds to System B, 6 = Aofb corresponds to System C, and g{t) = g{92{t)). 
Notice that each trajectory reaches after a time r = 2h/uj2. 

Motivated by the observation that z{6i) = 0{b^) in the relevant rectangle, we rescale the 
variable 6hy 6 = b'^6, and define z{6i) = ^z{6i). This gives 

ei = -uj,-{b^as + 6)z{ei)g{t) (10) 
6=0. (11) 

To find a, we solve ([IOl)-([II]) over the time interval [0, r], starting from (6'i(0), 5(0)) = (1 - 6, 0); 
to find a', we do the same, starting from (1 — 6, fe^Aofb). Applying Gronwall's Lemma gives 
\a' — a\ < fo^Aofb exp(Lr), where L is the Lipschitz constant for the allied vector field. To 
estimate L, note that \g\ = C(i), \z\ = 0(1), and \z'\ = C(i). This gives L = 0{l/b), and 
exp(Lr) = 0{1). Therefore, \a' — a| < fc2fe^Aafb for some constant Note that this constant 
can be made independent of tOi,uj2, as or Ofb- 

Recall from the first part of the proof that to obtain the desired result for System C, it suffices 
to guarantee |a' — a| < h. This happens when 

Aafb < 7^ := Kb~'Au . (12) 

□ 

Lemma l33t There exist b2, C > andxi G (0, 1—6) such that for all admissible (cui, ci;2, ag, Ofb), 
ifb < 62 cind Aflfb > CAcj, then T{xi) > xi. 

Proof. All orbit segments considered in this proof run from fl {^^i G (0, 1)} to We 
assume afj > 0; the case afj < is similar. First, we fix xq, si > b so that for all admissible 
(c<;i,c<;2, Off, ttfb), the trajectory starting from xi intersects H = {1 — 6 < 62 < 1 + b} in H (1 
{01 G (1 + Xq, |)}. Such xq, Xi clearly exist for small enough b. Starting from xi, we compare 
the trajectories for Systems A, B and C. We know that the trajectory for System A will end in 
(1 + Xi, 1 + 6). Thus to prove the lemma, we need to show the trajectory for System C ends to the 
right of this point. This comparison is carried out in two steps: 

Step 1: Comparing Systems A and B. We claim that the horizontal separation of the end points 
of these two trajectories is < cAu for some constant c > 0. It is not a necessary assumption, but 
the comparison is simpler if we assume b is small: First, a separation ciAco in the 62 direction 
develops between the trajectories as they flow linearly from (xi, 6) to = (1 — 6). Next, while 
the trajectories flow through V, Gronwall's lemma can be used in a manner similar to the above 
to show they emerge from V with a separation < C2AUJ in the 62 direction. Third is the region 
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of linear flow, resulting in a separation < caAo; in the 9i direction as the trajectories enter H. 
Gronwall's lemma's is again used in the final stretch as the trajectories traverse H. 

Step 2: Comparing Systems B and C. Notice that up until they reach Ei_6, the two trajectories are 
identical. In H, their 62 coordinates are equal, and the crossing time is r = — . Let Of{t) and 
Of{t), t e [0, r], denote their 9i coordinates while in H. We write 

0?{r)-ef{T) = r aii{z{e'^{t))-z{9f{t)M02{t))dt + r Aaii,z{9^{t))g{e2{t))dt. 
Jo Jo 

The first integral is > by design: via our choice of xi, we have arranged to have 1 < 9i{t) < 
(^) < 1' the 2;-function is monotonically increasing between 9i = 1 and 6'i = |. As for 
the second integral, we know z{9f{t)) is bounded away from since 1 + xq < 9f{t) < |, so the 
integral is > dAafb for some constant d > 0. It follows that T{xi) > Xi if dAofb > cAo;. □ 
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